clc

tic
format short

%% Set exogenous parameters
global eS eB gammaR1 gammaB A3g omega3g omega3b delta xihigh xilow cf1 cf2 cf3 ceq  wE wB wS ...
       rho  mc cI phiI CC phiD cud cud2 CR LCR LevR NSFR LR LIQreg scale Rout wN...
       InCR InLCR InLevR InNSFR InLR InLIQreg InScale IntD IntLIQ IntI IntE Dscale LIQscale Iscale Escale InLIQS

%% Parameters
% % 
eS= 2.5; %endowment of savers
eB= 0.1; %endowment of banker
cf1=1; %time discount savers t=1
cf2=0.7; %time discount savers t=2
cf3=cf2^2; %time discount savers t=3
gammaR1=1; %risk aversion savers
gammaB=1; %risk aversion banker
rho=0.15; %coefficient multiplying t=1 utility of banker (gamma in the paper---reserved gamma for risk aversion in code)
omega3g=0.8; %probability of good state
omega3b=1-omega3g; %probability of bad state
xilow=0.01; %lower value for xi
xihigh=2.2; %higher value for xi
delta=0.5; %portion of impatient savers
cI=0.1; %coefficient multipying effort cost
phiI=3; %coefficient in effort cost function
mc=0.1; %monitoring cost
ceq=0.27; %coefficient multipling payments services from deposits
phiD=1; %coeeficient in payment services function
cud=1; % auxilary variable to change payment services function
cud2=0; % auxilary variable to change payment services function
A3g=4.5; %productivity in good state
Rout=1; %Outside option of the banker

% Welfare weights
wE=0.0; % weight on entrepreneur---zero gives PE---adjust to obtain planner's solution
wS=0.0; % weight on depositor---zero gives PE---adjust to obtain planner's solution
wB=1; %weight on banker---keep at one for both PE and SP

%Other
wN=1; %weight on deposits in NSFR
CC=1; %Index is equal to one if deposit contract is complete; otherwise zero
InLIQS=1; %Index whether in equilibrium LIQ_1^S=0---set to one if CE(13)<0 in code

%% Indices for active regulation (for Table 2 set all to zero)
       % Initially set all regulations inactive (zero value for indices)
       InCR=0; %capital requirement
       InLCR=0; %liquidity coverate ratio
       InLevR=0; %leverage ratio
       InNSFR=0; %net stable funding ratio
       InLR=0; %liquid to loans ratio
       InLIQreg=0; %liquid to assets ratio
      % Also set indices whether Pigouvian regulation is active to zero 
       InScale=0; %Pigouvian tax to target I-E
       IntI=0;%Pigouvia tax to target I
       IntD=0; % Pigouvian tax to target D
       IntE=0; %Pigouvian tax to target E
       IntLIQ=0; %Pigouvia tax to target LIQ

%% Define global variables that will be used in later files
 scale=0;
 Dscale=0;
 LIQscale=0;
 Iscale=0;
 Escale=0;
 CR=0;
 LevR=0;
 LCR=0;
 NSFR=0;
 LR=0;
 LIQreg=0;

%% Initial values for Private Equilibrium
CEinit=[0.862211969545343;0.0516328946182285;0.875460152643298;0.0383847115202732;0.900670624961223;3.09726477083171;0.717302758273125;-9.34514822129717e-28;2.43446072721634;0.277308381893769;0.895309382147000;1.31951201549580;-0.0561324215508608;0.819705661331778;-1.03454495457085e-58;0;0;0;0;0;0;0;0;0;0;0];

%% Implementation of planner's solution in Tables 5 
    % % Planner's solution for wE=0 and wS=0.2
    % %Cand_Table5=[I in SP; LIQ in SP; D in SP; E in SP]
      Cand_Table5=[0.785488634413642;0.220554382366427;0.961977767159920;0.0440652496201495];
    % % Use Cand_Table5 to compute ratios that regulation targets
                CRcandT5=Cand_Table5(4)/Cand_Table5(1);
                LIQregcandT5=Cand_Table5(2)/(Cand_Table5(2)+Cand_Table5(1));
                DscalecandT5=Cand_Table5(3);

%% Create vectors for different regulations in Table 5
InCRT5=[0,1,1,1];
CRT5=[0,CRcandT5,0.0567,CRcandT5]; %adjusted CR so that CR is bidning by iteslf or both regulations are binding when CR and ell are used jointly but without Pigouvian tax
InLIQregT5=[0,0,1,1];
LIQregT5=[0,0,0.2098,LIQregcandT5]; %adjusted LIQreg (ell) so that both regulations are binding when CR and ell are used jointly but without Pigouvian tax
IntDT5=[0,0,0,1];
DscaleT5=[0,0,0,DscalecandT5];

%% Create vector to store utilities
UB=zeros(1,4);
USstar=zeros(1,4);
UEstar=zeros(1,4);

%% Create matrix to store PE and SP equilibrium variables reported in Table 2
EqSol=zeros(17,4);

%% Solve for equilibrium---looping over wST2 and wET2
options=optimset('Display','iter','Tolfun',1.0e-12, 'MaxFunEvals',100000,'MaxIter',10000);

for i=1:size(CRT5,2)


    InCR=InCRT5(i);
    CR=CRT5(i);
    InLIQreg=InLIQregT5(i);
    LIQreg=LIQregT5(i);
    IntD=IntDT5(i);
    Dscale=DscaleT5(i);

[CE,F]= fsolve (@code_final_submission,CEinit,options);

%% Equilibrium Values
I=CE(1); %loan to entrepreneur
L=CE(2); %liquid asset holdings at t=1
D=CE(3); %deposits
E=CE(4); %equity
xistar=CE(5); %run threshold
rI=CE(6); %loan rate (net). Gross rate in the paper is R_I=1+rI
r3D=CE(7); %late deposit interest rate (net). Gross rate in the paper is R_D=1+r3D
r2D=CE(8); %early deposit interest rate (net). Gross rate in the paper is r_D=1+r2D
psiBS=CE(9); %L multiplier on balance sheet constraint
psiGG=CE(10); %L multiplier on global game constraint
psiLD=CE(11); %L multiplier on loan demand schedule
psiDS=CE(12); %L multiplier on deposit supply schedule
LIQ1R=max(CE(13),0); %Savers' liquid asset holdings
slackSP=CE(14)^2; %Slack participation constraint of banker--for planner's problem
psiSP=CE(15); %L multiplier on participation constraint of banker
lambdaCR=CE(16); %L multiplier on capital requirement
lambdaLCR=CE(17); %L multiplier on liquidity coverage ratio requirement
lambdaLevR=CE(18); %L multiplier on leverage ratio requirement
lambdaNSFR=CE(19); %L multiplier on net stable funding ratio requirement
lambdaLR=CE(20); %L multiplier on liquidity to loans ratio requirement
lambdaLIQreg=CE(21); %L multiplier on liquidity to assets ratio requirement
tscale=CE(22); %Pigouvian tax on scale implemented as a tax on deposit interest payment
tD=CE(23); %Pigouvian tax on deposits
tLIQ=CE(24); %Pigouvian tax on liquid asset holdings
tI=CE(25); %Pigouvian tax on loans
tE=CE(26); %Pigouvian tax on equity

%% Out of equilibrium withdrawals
thetastar=(L+xistar*I)/(D*(1+r2D));
hatl=((xistar*I+L)*(1+rI)-xistar*(D*(1+r3D)+mc/omega3g))/(D*((1+r2D)*(1+rI)-xistar*(1+r3D)));
ls=[thetastar;hatl;delta];

%% Probability of a run
q=(xistar-xilow)/(xihigh-xilow);

%% Regulatory ratios
CR=E/I;
LevR=E/(I+L);
LCR=(L+xilow*I)/(D*(1+r2D));
NSFR=(E+wN*(1-delta)*D)/I;
LR=L/I;
LIQreg=L/(L+I);

%% Preference for liquidity parameteters
if phiD==1
gD=cud2+(log(cud+D*(1+r2D)))*ceq;
else
    gD=cud2+((cud+D*(1+r2D))^(1-phiD)/(1-phiD))*ceq;
end
dgD=(1+r2D)*(cud+D*(1+r2D))^(-phiD)*ceq;
d2gD=-phiD*(1+r2D)^2*(cud+D*(1+r2D))^(-phiD-1)*ceq;

%% Utilities
UB(i)=rho*U1(eB-E,gammaB,cf1)+integral(@(xi)omega3g*((xi*I-delta*D*(1+r2D)+L)/xi*(1+rI)-(1-delta)*D*(1+r3D))-mc,xistar,xihigh,'ArrayValued',true)/(xihigh-xilow);
UEstar(i)=(xihigh-xistar)/(xihigh-xilow)*(cI*(phiI-1)*I^phiI);
USstar(i)=U1(eS-D-LIQ1R,gammaR1,cf1)+MU1(eS-D-LIQ1R,gammaR1,cf1)*(D+LIQ1R)+(xihigh-xistar)/(xihigh-xilow)*(gD-dgD*D);

%% Expected dividends
EDIV=integral(@(xi)omega3g*((xi*I-delta*D*(1+r2D)+L)/xi*(1+rI)-(1-delta)*D*(1+r3D))-mc,xistar,xihigh,'ArrayValued',true)/(xihigh-xilow);

%% Check that early deposit rate r2D is zero: Show that first-order condition of banker with respect to r2D evaluated at the equilibrium values is negative

UBdr2D=-omega3g*delta*(1+rI)*D*log(xihigh/xistar)/(xihigh-xilow);
hatldr2D=-hatl*(1+rI)/(((1+r2D)*(1+rI)-xistar*(1+r3D)));    
GGdr2D=hatldr2D*omega3g*D*(1+r3D)-D*(thetastar-delta);    
LDdr2D=-omega3g*(A3g-(1+rI))*(delta*D/I*log(xihigh/xistar)/(xihigh-xilow));
d2gDr2D=-phiD*D*(1+r2D)*(cud+D*(1+r2D))^(-phiD-1)*ceq;
DSdr2D=(1-q)*cf2*delta+(1-q)*d2gDr2D+(1-q)*(cud+D*(1+r2D))^(-phiD)*ceq;     
 % % First order condition
FOCr2D=UBdr2D+psiGG*GGdr2D+psiLD*LDdr2D+psiDS*DSdr2D;

%% store equilibrium variables for PE and SP equilibria
EqSol(:,i)=[I;L;D;E;1+rI;1+r3D;q;LIQreg;LevR;1+rI-(1-delta)*(1+r3D);I+L;I-E;EDIV;UEstar(i)-UEstar(1);...
        USstar(i)-USstar(1);UB(i)-UB(1);(UB(i)+USstar(i)+UEstar(i))-(UB(1)+USstar(1)+UEstar(1))];

end

%% Create string with variables names
variables={'I';'L';'D';'E';'R_I';'R_D';'q';'\ell';'k';'R_I-(1-\delta)R_D';'I+L';'I-E';'E(Div)';'\Delta U_E';'\Delta U_S';'\Delta U_B';'\Delta S_{sp}'};

%% Produce Table 5
Table5=table(string(variables),EqSol(:,1),EqSol(:,2),EqSol(:,3),EqSol(:,4),'Variablenames',{'Var','PE','CR','CR & ell','CR,ell &tau_D'});

%% Computing time
toc

t=toc;